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Abstract 

A one-dimensional version of the second-order transition model based on the sheared flow amplifi- 
cation by Reynolds stress and turbulence supression by shearing is presented. The model discussed 
in this paper includes a form of the Reynolds stress which explicitly conserves momentum. A 
linear stability analysis of the critical point is performed. Then, it is shown that the dynamics of 
weakly unstable states is determined by a reduced equation for the shear flow. In the case in which 
the flow damping term is diffusive, the stationary solutions are those of the real Ginzburg-Landau 
equation. 
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I. INTRODUCTION 



The existence of a shear flow layer at the tokamak edge and in ohmic discharges has been 
known for a long time . It was later found that similar edge shear flow layers existed in 
other confinement devices. This seems to be a generic feature of confined plasmas. In the 
last years, a great deal of attention has been directed to the formation of shear flow layers 
and the corresponding region of radial electric field gradient to understand the improved 
confinement regimes. Many of the recent theoretical developments in this direction have 
been focused on barrier formation 3] or zonal flows 4-J. On the experimental side, much 
progress has been done in the visualization of edge turbulence and flows ([^j, including 
specific applications to the emergence of the shear flow layer in the TJ-II stellarator (3)), 
in which we are especially interested. 

Here, we want to turn back to the basic plasma edge shear flow layer. In stellarators, 
unlike in tokamaks, one can operate at densities for which no shear layer is present in the 
plasma edge, being thus possible to study its formation. 

The emergence of theplasma edge shear flow layer as the density increases in the TJ-II 



stellarator 



is shown 



transition. It is consistent 
fication by turbulence 1121 
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ent iii 

3,13 



to have the characteristic properties of a second order phase 
111 with a simple transition model that couples shear flow ampli- 

n 

with turbulence suppression by sheared flows [14|. The model 
used in interpreting the TJ-II results is based on a transition model [l5] initially introduced 
to explain the transition from the low confinement mode (L mode) to the high confinement 

n 

mode (H mode) [16( in magnetically confined plasmas. This model consists of two envelope 
equations for the fluctuation level and mean poloidal flow. A later extension of the model 



171 ] included a third equation to account for the pressure gradient contribution to the radial 



electric field. This second model shows the existence of two critical points, the second one 
causing the first order transition that has been associated with the L to H transition. For 
this transition there is a hysteresis cycle and the transition is characterized by an S-curve 
(see for instance 18]). The first critical point leads to a second order transition (conse- 
quently, it does not have a hysteresis cycle), which has been identified with the emergence 
of the plasma edge shear flow layer. 

In this paper, we focus on this second order transition. By excluding the diamagnetic 
term in the momentum balance equation, only this transition is included in the model. This 
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simplification is reasonable because there is a large range of densities separating the two 
critical points. For the same reason, the range of plasma parameters considered is not yet 



in an L-mode regime. Therefore, we do not expect avalanche-like transport 13, |20[ and 



the transport terms can be represented by purely diffusive terms. Concretely, we discuss an 



extension to 1-D of the original model used in comparison with t 



re experimental data 



21 



221 ] . here we formulate 



In contrast with previous 1-D extensions of the transition model 
the Reynolds stress term as a momentum conserving term. The model is defined by three 
1-D partial differential equations describing the evolution of the turbulent fluctuation level 
E, the averaged poloidal velocity shear U and (minus) the pressure gradient N. It predicts a 
second-order phase transition with order parameter U and control parameter T, the particle 
flux, which enters the model through the boundary conditions. For T below a critical value 
r c the stable stationary solutions have U — 0, whereas for T > T c such solutions are unstable 
and undergo a transition to states with U ^ and reduced turbulence fluctuations. 

After performing a detailed stability analysis of the model we study two interesting special 
cases depending on the form of the flow-damping term: collisional drag and collisional 
diffusion. For both of them we find reduced equations describing the dynamics of weakly 
unstable states. In the latter case, the reduced equation is closely related to the Ginzburg- 
Landau equation for second-order phase transitions. In particular, the stationary solutions 
of our equation are exactly those of the Ginzburg-Landau one. 

The paper is organized as follows: 

In Section II we introduce the one-dimensional transition model with momentum con- 
servation. Section III is devoted to the study of the fixed points of the model, their linear 
stability properties and a general discussion of the critical conditions. In Section IV we 
consider the dynamics near the critical point. Section V contains the conclusions and an 
outline of future research lines. 



II. THE ONE-DIMENSIONAL TRANSITION MODEL 

The relevant plasma edge region to which this model is applied corresponds tor G [r , a], 
(a — ro) /a ~ 0.1, where a is the minor radius of the plasma. We take the slab geometry 
approximation in representing this region and use as coordinate x := (r — r )/(a — r ), so 
that x G [0, 1]. 
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The fields of our model will be the fluctuation level envelope £ := ((nk/no) 2 ) 1 / 2 , the 
averaged poloidal shear flow IA := d(Vg)/dr and (minus) the averaged pressure gradient 
M := —d(p)/dr, where (•) denotes ensemble average. A suitable one- dimensional general- 



ization of the model discussed in Ref. ll| requires a form of the Reynolds stress which con- 
serves momentum. Using a pressure-gradient-driven turbulence model and assuming densely 
packed turbulence, a quasi-linear calculation yields the following form for the Reynolds 
stress: 

(V x Ve) = a 3 £ 2 d x (V e ) + D 2 £ 2 d 3 x (V e ), (1) 



which is very similar to the expression previously derived in Ref. |23j in the context of zonal 
flow dynamics. Notice that measures the strength of the Reynolds stress and non-zero 
D 2 is needed for the spectrum of the instability to be bounded. 
The model discussed in the present work is: 

dtS = jJ\fS - a x £ 2 - a 2 U 2 £ + d x [(D + D 1 £)d x 8] (2a) 
dtU = -fitU + fi 2 d 2 x U - a 3 d 2 x {£ 2 U) - D 2 d 2 x {£ 2 d 2 x U) (2b) 
d- t N = d 2 x [(D 3 £ + ^4)^] . (2c) 

Here, 70AA is the linear growth rate of the characteristic instability and a\ is computed 
from the nonlinear saturation condition of the instability, both in the absence of sheared flow. 
The a 2 coefficient is derived from the condition of turbulence supression by sheared flow. 
Do an d -D 4 are neoclassical diffusivity coefficients, whereas Di and D 3 are the coefficients of 
anomalous diffusivity multiplying the fluctuation level. Finally, Ji\ and Jx 2 are the coefficients 
of the collisional flow-damping terms. 

We can eliminate the explicit dependence on the parameters 70, «i and a 2 by means of 
the following change of variables: 

E := —£, U := A l^-U, t := j i, N := J\f, (3) 
7o V 7o 

and Eqs. ([2]) read in terms of E(x,t), U(x,t) and N(x,t): 

d t E = NE — E 2 — U 2 E + d x [(D + DiE)d x E] , (4a) 
d t U = -mU + fi 2 d 2 x U - a 3 d 2 x (E 2 U) - D 2 d 2 x (E 2 d 2 x U), (4b) 
d t N = d 2 x [(D 3 E + ,D 4 )iV] , (4c) 
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with D = A)/7o, Di = Di/ax, D 3 = D 3 /ax, D A = D^/jo, fix = /ix/jo, A*2 = ^2/70, 
a 3 = a^o/af, D 2 = -D270/ «i- Finally, we choose the boundary conditions: 

d x U(u, t) = d 3 x U(u, t) = 0, u = 0, 1, Vt, (5a) 

d x E(u,t) = 0, u = 0,1, Vt, (5b) 

(L> 3 £ + ^4)^V| (0)f) = T, d x N\ {lt) = 0, (5c) 

where V is the particle flux and is the natural control parameter of the model. 



III. FIXED POINTS AND LINEAR STABILITY ANALYSIS 



It is obvious from Eq. ( 14bl) that any fixed point of the model must satisfy U — 0. Then, 
from Eq. (T4"a|) we find that NE = E 2 . It only remains to use the boundary condition of N 
at x = and we finally have that 

(i) There always exists a fixed point 

U f = 0, E f (T) = Nf(T) = J^- D 4 2 + 4D 3 r - £4) • (6) 

(ii) If D4 7^ there exists a second fixed point (always unstable, see below) 

u' f = 0, E' f = 0, iv}(r) = t/d 4 . (7) 

Let (So, U , N ) be a fixed point and linearize the equations (j4j) around it: 

E(x,t) =E + ^ t+lkx , 
U(x, t) = Cue" t+ikx , 
N(x,t) = N + i N e< t+tkx . 

For the fixed point (E'f, U'f, N'f) the eigenvalue condition leads to the dispersion relations: 

7 - N' f + D k 2 = 0, 
7 + fix + fi 2 k 2 = 0, 

7 + L^ 2 = 0. (8) 

It is clear that this fixed point is unstable for any T > 0, since the first dispersion relation 
gives, for k = 0, 7 = T/D 4 . 



The fixed point (Ef, Uf, Nf) is more interesting. In this case we have: 
7 = -n x - fi 2 k 2 + a 3 E 2 k 2 - D 2 E 2 f k 4 , 



J f ~^2^ ft 

(7 + E f + (D + 2D 1 E f )k 2 )( 1 + (D 4 + D 3 E f )k 2 ) + D 3 E 2 f k 2 = 0. (9) 

The second equation does not give solutions with 7 > 0. However, the first one can yield an 
instability. Then, we ask under which conditions 

j(k) = -/n - fi 2 k 2 + a 3 E 2 f k 2 - D 2 E 2 f k* (10) 

is positive. The neutral modes, i.e. the values of k for which ^(k) = are given by: 



a 3 E 2 -fi 2 ± J(a 3 E 2 - fi 2 ) 2 - 4^D 2 E 2 

kl = ¥ (11) 

± 2D 2 E 2 ' V ; 

which has real solutions k-, k + if and only if 

a 3 E 2 -fi 2 - 2E f ^JT 2 > 0. (12) 
Now observe that the boundary conditions imply the quantization of k. Namely, 

k = nvr, < n G Z. (13) 



Therefore ( 1121) is only a necessary condition for the existence of instabilities. In addition, 
there must exist some k = nir such that fc_ < nir < k + . The critical point is defined by the 
minimum value of the flux, r c , for which there exists an unstable mode k c = n c ir. 
We define E c := E f (T c ). From ([10]) we find that 

2 = /X! + /i 2 fc c 2 

c kl{a 3 - D 2 k 2 ) ■ 1 1 

Thus, in particular, if a 3 /D 2 < it 2 there are no unstable modes, no matter how much we 



increase T. If a 3 /D 2 > n 2 at least k = n can become unstable. Actually, if y / a 3 /D 2 G 
(it, nn), there exist n — 1 potentially unstable modes k = it, 2tt, . . . , (n — l)ir. 

IV. DYNAMICS NEAR MARGINAL STABILITY 

Our aim is to find approximate equations for the dynamics of Eqs. (j4j) near (and above) 
the critical point. To that end we perform an expansion with parameter 



£-1 (15) 

-L r. 



for small 5. Explicitly, we take: 



E = E C + 5 2 E 2 + ... (16a) 

U = 5U! + ... (16b) 

N = N C + 5 2 N 2 + ... (16c) 

r = r c (i + «5 2 ) (i6d) 

(recall that N c = E c ) and we perform a rescaling of the coordinates: 

fj = Sx, f = SH. (17) 
Expanding the equation of N gives: 

N 2 = E c (l- D " E2 ) (18) 

c \ d a + d 3 eJ y } 



whereas the equation for E yields: 



N 2 — E 2 — Uf = (19) 



and using ([TBI) we obtain for E 2 : 



E, = Di + DiEc (E c - Ul) . (20) 
2 L> 4 + 2D 3 E C K 1 ; V ; 

Consequently, for weakly unstable states, the problem of studying the dynamics of our 
model consists in finding an approximate equation for the dynamics of U, E and N being 
determined at the end of the day from the slaving conditions (1201) . ([TBI . The form of the 
reduced equation for U is quite different for \i\ ^ and \i\ = 0. That is why in the next 
sections we study the cases fi 2 = (collisional drag) and fix = (collisional diffusion) 
separately. If the collisionality at the plasma edge is high enough, the damping is essentially 
diffusive and the limit \i\ = is a good approximation. At low collisionality the magnetic 
pumping dominates and the relevant limit is /j, 2 = 0. 

A. Collisional drag 

Let us set \i 2 = 0, jii ^ 0. The dispersion relation reads: 

1 (k) = -fi 1 + E 2 f k\a 3 -D 2 k 2 ) (21) 
7 



so that 

Ec := yjw^-D 2 ki) ■ (22) 

As pointed out above, if a 3 /D 2 < it 2 there are no unstable modes. If it 2 < a 3 /D 2 < 47r 2 
only k = tt can be unstable and is, of course, the critical mode. If a 3 /D 2 > 4-7T 2 there are at 
least two modes which may become unstable and which of them is the critical mode depends 
on the quotient a 3 /D 2 . Hence, in general, k = n is not the most unstable mode (see Figs. 
Hand 125. 

In order to find a reduced equation for the weakly non-linear dynamics of U we expand 
(14bp keeping terms up to order 5 7 : 

where we have made use of (120|) . Finding a reduced equation without an explicit dependence 
in 5 seems difficult in this case. However, we can still obtain a useful reduced equation by 
simply taking ()23l) and going back to the original variable U and coordinates x, t: 

d t U = -^17 - E f {Vf (a 3 dlU - D 2 dtU) + 2 ^+^ [a 3 d 2 x U 3 + D 2 d 2 x (U 2 d 2 x U)} . 

(24) 

The time-evolution predicted by this equation is compared to the original model in Fig. 
E] for 5 = 0.274. The dynamics is very sensitive to 5 and Eq. (|24[) ceases to describe it 
accurately for larger values of the expansion parameter. However, it is remarkable that 
even for values of S of order 1, the stationary solutions of the reduced equation give good 
approximations of the exact ones (see Fig. HI where 5 = 0.82). 

B. Collisional diffusion 

In this section we take n\ = and fi 2 7^ 0. In this case 

7 (fc) = (-0a + a 3 E 2 f - D 2 E)k 2 )k 2 (25) 
and the critical point is given by 



As we already know, if a 3 / D 2 < it 2 there are no unstable modes. Unlike the case of collisional 
drag, if a 3 / D 2 > 7r 2 the most unstable (i.e. critical) mode is always k = it (see Fig. [5]). In 
addition, for typical values of a 3 and D 2 , a 3 ^> ir 2 D 2 , and we can use the approximation 

E C = M. (27) 

V "3 

We are now ready to derive a reduced equation for the dynamics of weakly unstable states 
in this case. The expansion of (14bl) in powers of 5 yields: 



8 b d f U x = 5 3 d 2 fj [(/x 2 - a 3 E 2 c )Ut - 25 2 a 3 E c E 2 U l - 8 2 D 2 E 2 d 2 U 1 ] . (28) 

Noting that the first term on the right-hand side vanishes if we use the approximation (12711 . 
we are left with 

d f U x = -d 2 [2a 3 E c E 2 U 1 + D 2 E 2 c dp x ] . (29) 
Finally, we use (120]) and perform one more change of variables: 



° = ^ Ul > T = ^{ D i + 2D 3 E C ) T > ^=^D 2 (D 4 + 2D 3E y (30) 



^E, A 



obtaining the definitive form of the equation describing the weakly non-linear dynamics of 
U: 

d T a = -d 2 [a-a 3 + d 2 a}. (31) 

In Fig. [6] we show a comparison between the dynamics of Eq. (|3~1|) and that of the original 
model, whereas Figs. [7] and [8] are comparisons of the stationary solutions. 

A remark is in order at this point. A consequence of using the approximation ( 1271) in 
Eq. (128j) is that the linear growth rate is modified. An easy calculation shows that Eq. (13~TT) 
gives the actual linear growth rate (at order 5 2 ) only if we make the replacement 

5 2 » S 2 + ^ (32) 

«3 

which is natural, since the accuracy of our approximation depends on the quotient D 2 ti 2 /a 3 . 
Thus, very close to the critical point one has to take into account the above correction of 
S when comparing the results of ( 131]) with those coming from the integration of the exact 
equations of the model. 

Using that at the boundary, d x a = d 3 a = 0, and taking an initial condition such that 
Jq U(x)dx = 0, the stationary solutions of Eq. (13T1) are exactly the solutions of: 



a — a 



3 1 a2 



dla = 0, (33) 



which is the time- independent Ginzburg-Landau equation for second-order phase transitions. 
It is the same equation as in the case of non- momentum-conserving Reynolds stress with 



a collisional drag (see Ref. (24|). We follow the lines of Ref. [24j to obtain analytical 
expressions for the stationary solutions. The key observation is that there exists a 'conserved 
quantity'. Namely, 

^(<V) 2 + y(i-y) = C, CeR. (34) 
This allows to reduce the solutions to quadratures: 

f da f (35) 

J vV 2 -6 + )(a 2 -M J V2 

with b± — 1 ± y/l — AC. Defining m = b^/b + and performing the change of variable 
a = y^Tsin0 we can recast the solutions in the following form: 

T^= / / d9 .2n ~ A > (36) 

z Jo yl — m sin 9 
where A is an integration constant. 

Hence, the solution for a(rj) may be expressed as 

a(rj) = ^/blsn (^]j^V + (37) 

where sn(y\m) stands for the Jacobi elliptic function, which is periodic in y. Its period is 
P = AK(m), with 

-tt/4 dQ 




K(m) = 4 / (38) 
Jo yl - msin 9 

the complete elliptic integral of the first kind. In order to write explicitly the solution in 
terms of the original variable U and coordinate x, recall that 



U = 5VE ef r, V = S s l D2{Di + 2D3Ec) x. (39) 



Then, 



(40) 



Thus, the wavelength of the solutions is 



AKim) I D 2 (D 4 + 2D 3 E C ) 

* Xl b + a 3 (D, + D 3 E c )- 141 ] 
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The boundary conditions determine the integration constants C and A. In particular, 
they imply the quantization condition 

n^ = l, 0<neZ. (42) 

V. CONCLUSIONS AND FURTHER WORK 

We have introduced a one-dimensional version of the phase transition model considered 
in Ref. 



ll| including a Reynolds stress term with manifest momentum conservation. The 
model consists of three envelope equations for the fluctuation level E, the poloidal shear 
flow U and the density gradient N. It possesses a critical point corresponding to a second 
order transition whose natural control parameter is T, the particle flux. Below the critical 
value, T c , the model has a non-trivial fixed point with U = 0. Through a linear stability 
analysis we have shown that if T > T c the stationary solutions have non-zero shear flow and 
reduced turbulent fluctuations. We have also studied the dynamics of the model near (and 
above) the critical point. Defining a suitable expansion around the critical point we have 
derived slaving conditions for E and N, so that they are determined from the value of U. 
Then, we have found reduced equations for the weakly non-linear dynamics of U. In the 
case of diffusive shear flow damping the reduced equation is related to the Ginzburg-Landau 
equation, which allows us to work out analytical expressions for the stationay solutions of 
weakly unstable states. 

The analysis performed in this work shows some interesting differences with respect to 
previous one-dimensional versions of the model in which momentum conservation is not 



implemented (see Ref. 24j), the most relevant of them concerning the nature of the fixed 



points and the instabilities. 



In Ref. 24| the most unstable mode is always k — 0, which is related to the fact that in 
the model studied therein there exist two (non-trivial) fixed points, one of them with U = 
and the other one with U ^ 0. The introduction of a momentum-conserving Reynolds stress 
makes the fixed point with non-zero U disappear, so that all stationary solutions with non- 
zero shear flow have non-trivial spatial structure. This is connected to the results derived 
in previous sections regarding the linear stability analysis of the critical point. In the model 
with momentum conservation k = is always stable and the discussion of the structure of 
the most unstable (critical) mode is more complicated. In general, the critical mode depends 

11 



on the values of the parameters. The presence of non-zero critical mode is in agreement with 
the experimental findings reported in Ref. [lo| . 

A next step in the development of the model would be to incorporate the diamagnetic term 
in order to study the L to H transition. In the L-mode regime the model must incorporate the 
transport mesoscale, which can be achieved by formulating transport equations in terms of 
fractional derivative operators 



of a react ion- diffusion system 



] . Such modification of the model will lead to the dynamics 
26|. These issues will be addressed in future publications. 
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Figure 1: 7 as a function of T near the critical point. The values of the parameters are 03 = 10~ 2 , 
Z?2 = 10 -4 , D3 = 10~ 2 , Dq = D\ = D4 = 0, n\ = 1, fj,2 = 0. The critical mode is k = 2ir and 
the critical point is given by r c = 0.041853. k = 3ir becomes unstable at T = 0.1008 and k = tt at 
r = 0.1124. 
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Figure 2: j(k) as a function of k near the critical point. The values of the parameters are = 10~ 2 , 
L>2 = 10~ 4 , D3 = 10~ 2 , Do = D\ = Z?4 = 0, /ii = 1, fj,2 = 0. The critical mode is A; = 2tt and the 
critical point is given by r c = 0.041853. 
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Figure 3: Time-evolution of U(0) computed from the integration of the original model (solid) and 

from the reduced equation (I24j) (dashed). The values of the parameters are 03 = 

D 3 = 10" 2 , D = Di = D A = 0, m = 1, fiQ = 0, r = 0.045. The critical point is T c = 0.041853. 
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Figure 4: Stationary solutions computed from the integration of the original model (solid) and 
from the reduced equation (I24p (dashed). The values of the parameters are 03 = 1(T 2 , D 2 = 10~ 4 , 
D 3 = 1(T 2 , D Q = D x = D 4 = 0, m = 1, iii = 0, r = 0.07. The critical point is T c = 0.041853. 
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Figure 5: "f(k) as a function of k near the critical point. The values of the parameters are 03 = 10 
D 2 = 1CT 4 , D 3 = 1CT 2 , D = Di = D A = 0, m = 0, \i 2 = 1- The critical point is T c = 1.1095. 
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Figure 6: Time-evolution of U(0) computed from Eq. (|31|) (dashed) and from the integration of 
the original model (solid). The values of the parameters are 03 = 10~ 2 , D2 = 10~ 4 , D3 = 10~ 2 , 
Dq = D\ = D4 = 0, jjL2 = 1, jJL% = 0, T = 1.2. The critical point is T c = 1.1095. 
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Figure 7: Structure of the bifurcation at the critical point. The values of the parameters are 
a 3 = 10~ 2 , D 2 = 10" 4 , D 3 = 10~ 2 , D = D x = D 4 = 0, // 2 = 1, m = 0. The critical point is 
T c = 1.1095. The solid curve corresponds to the analytical solutions of Eq. [33l The points have 
been computed from the numerical integration of the original model. As we can see, the reduced 
equation gives very accurate results even at values of 5 of order one. 
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Figure 8: Stationary solutions computed from the integration of the original model (solid) and 
from the GL equation (|33|) (dashed). The values of the parameters are 03 = 1(T 2 , D 2 = 10" 4 , 
D 3 = 1(T 2 , D = Di = D 4 = 0, ix 2 = 1, Ml = 0, T = 1.5. The critical point is T c = 1.1095. 
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